#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*http://gmplib.org/*/
#include <gmp.h> 

/* Installation instructions for GMP are available here:
http://gmplib.org/manual/Installing-GMP.html#Installing-GMP*/

/*Important note: All programs using GMP must link against the libgmp library,
for example: gcc myprogram.c -o myprogram -lgmp*/	

static const int terms = 50;
int x = 1;
int work_unit = 30;


/* prototypes */
void get_pi();
void factorial(unsigned long int n, mpz_t fact);


void factorial(unsigned long int n, mpz_t fact)
{
	unsigned long int  num = n;
	mpz_fac_ui (fact, num);
}

void get_pi()
{	
	int limit = work_unit * terms;
	int n;
	
	unsigned long int two;
	unsigned long int n2;
	unsigned long int n2p;
	two = 2;
	
	mpf_t dividend_float, divisor_float, pi;
	mpf_init(dividend_float);
	mpf_init(divisor_float);
	mpf_init(pi);
	
	mpz_t dividend;
	mpz_init (dividend);
	
	mpz_t base;
	mpz_init (base);
	mpz_set_ui(base, two);
	
	mpz_t divisor_A;
	mpz_init (divisor_A);
	
	mpz_t divisor_B;
	mpz_init (divisor_B);

	mpz_t divisor_C;
	mpz_init (divisor_C);
	
	mpz_t divisor_T;
	mpz_init (divisor_T);
	
	mpz_t divisor_F;
	mpz_init (divisor_F);
	
	mpf_t partial_result;
	mpf_init (partial_result);
	
	mpf_t final_result;
	mpf_init (final_result);
	mpf_set_ui(final_result, 0);
	
	for(n = 0; n < limit; n++){
		n2 = 2*n;
		n2p = 2*n+1;
		
		factorial((2*n), dividend); 
		mpz_pow_ui(divisor_A, base, n2);
		factorial(n, divisor_B);
		mpz_pow_ui(divisor_B, divisor_B, two);
		mpz_set_ui(divisor_C, n2p);
		mpz_mul(divisor_T, divisor_A, divisor_B);
		mpz_mul(divisor_F, divisor_T, divisor_C);		
		mpf_set_z(dividend_float, dividend);
		mpf_set_z(divisor_float, divisor_F);
		mpf_div(partial_result, dividend_float, divisor_float);
		
		mpf_add(final_result, final_result, partial_result);
	}
	
	mpf_mul_ui (pi, final_result, two);
	gmp_printf ("Calculated value for Pi %.*Ff\n", 20, pi);
	
	/*clean space used by the objects*/
	mpz_clear(dividend);
	mpz_clear(base);
	mpz_clear(divisor_A);
	mpz_clear(divisor_B);
	mpz_clear(divisor_C);
	mpz_clear(divisor_T);
	mpz_clear(divisor_F);
	mpf_clear(partial_result);
	mpf_clear(final_result);
	mpf_clear(dividend_float);
	mpf_clear(divisor_float);
} 


int main()
{
  
  printf("*** Operation to determine the value of Pi ***\n");
  printf("-> Number of work units: %i\n", work_unit);
  printf("-> Number of iterations: %i\n", work_unit*terms);
  get_pi();
  return 0;
}

